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Electronic excitation energies of molecular systems from the Bethe-Salpeter equation: 
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We review the Bethe-Salpeter equation (BSE) approach to the calculation of electronic excitation 
energies of molecular systems. We recall the general Green's function many-theory formalism and 
give the working equations of the BSE approach within the static GW approximation with and 
without spin adaptation in an orbital basis. We apply the method to the pedagogical example of 
the H2 molecule in a minimal basis, testing the effects of the choice of the starting one-particle 
Green's function. Using the non-interacting one-particle Green's function leads to incorrect energy 
curves for the first singlet and triplet excited states in the dissociation limit. Starting from the exact 
one-particle Green's function leads to a qualitatively correct energy curve for the first singlet excited 
state, but still an incorrect energy curve for the triplet excited state. Using the exact one-particle 
Green's function in the BSE approach within the static GW approximation also leads to a number 
of additional excitations, all of them being spurious except for one which can be identified as a 
double excitation corresponding to the second singlet excited state. 
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I. INTRODUCTION 



Time-dependent density-functional theory 

(TDDFT) [l| within the linear response formal- 
ism [2tJ1 is nowadays the most widely used approach 
to the calculation of electronic excitation energies of 
molecules and solids. Applied within the adiabatic 
approximation and with the usual local or semilocal 
density functionals, TDDFT gives indeed in many cases 
excitation energies with reasonable accuracy and low 
computational cost. However, several serious limitations 
of these approximations are known, e.g. for molecules: 
too low charge-transfer excitation energies [5[, lack 
of double excitations [6[, and wrong behavior of the 
excited-state surface along a bond-breaking coordinate 
(see, e.g., Ref. 0). Several remedies to these problems are 
actively being explored, including: long-range corrected 
TDDFT [^ Q which improves charge-transfer excitation 
energies, dressed TDDFT @, [H [Tli which includes 
double excitations, and time-de pen dent density-matrix 
functional theory (TDDMFT) [l3-[ll| which tries to 
address all these problems. 

In the condensed-matter physics community, the 
Bethe-Salpeter equation (BSE) applied within the GW 
approximation (see, e.g., Refs-ilTHlQi) is often considered 
as the most successful approach to overcome the limi- 
tations of TDDFT. Although it has been often used to 
describe excitons (bound electron-hole pair) in periodic 
systems, it is also increasingly applied to calculations of 
excitation energies in finite molecular systems |2C)| - (3l| . In 
particular, the BSE approach is believed to give accurate 
charge-transfer excitation energies in molecules [2^, [3l[ , 
and when used with a frequency-dependent kernel it is in 



principle capable of describing double excitations J32ll33l |. 
In this work, we examine the merits of the BSE ap- 
proach for calculating excitation energies of the proto- 
type system of quantum chemistry, the II2 molecule. The 
paper is organized as follows. In Sec. |TT1 we give a re- 
view of Green's function many-body theory. In Sec. IIIH 
we give the general working equations for a BSE calcu- 
lation within the static GW approximation in a finite 
spin-orbital basis, and the corresponding spin-adapted 
expressions for closed-shell systems. In Sec. lIV| we apply 
the equations to the H2 molecule in a minimal basis and 
discuss the possibility of obtaining correct spin-singlet 
and spin-triplet excited-state energy curves as a function 
of the internuclear distance. Sec. FVl contains our conclu- 
sions. 



II. REVIEW OF GREEN'S FUNCTION 
MANY-BODY THEORY 

We start by giving a brief review of Green's func- 
tion many-body theory for calculating excitation ener- 
gies. For more details, see e.g. Refs. [13, ]^ [SJ. 



A. One-particle Green's function 

Let |iV) be the normalized ground-state wave function 
for a system of N electrons described by the Hamiltonian 
H. The time-ordered one-particle equilibrium Green's 
function is defined as 



iG{l,2) ^{N\f[i>{l)i'^2)]\N) 

=e(ti-t2)(iV|^(l)*^(2)|7V) 
-e{t2-h){N\^^{2)^{l)\N). 



(1) 
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Index 1 stands for space, spin and time coordinates 
(ri,cri,ii) = (xi,ii). T is the Wick time-ordering oper- 
ator which orders the operators with larger times on the 



left and is the Heaviside step function. The whole time- 
dependence is contained in *&(!) = e'^*i^(xi)e~*^*^ and 
*t(2) = e*^*2(|/t(x2)e~*^*^ the annihilation and cre- 
ation field operators in the Heisenberg representation. 

In the absence of external potential, the system is 
invariant under time translation, therefore the Green's 
function depends only on r = ti — ^2- By introducing the 
closure relation for excited states with A'' — 1 or A^ -I- 1 
particles, one can get 

iG(xi,x2;T) = 

0(r)^(iV|V3(xi)|iV + l,A)(iV + l,A|v3^(x2)|iV) 

A 

- ^(-^) E<^l^'^(^2)|iV - 1, 1){N - l,/|V;(xi)|iV) 



X e 



-i(EM-EM-i,i)'' 



(2) 



where E'^r, E'^r+i^^ and Epf-u are the energies of the 
ground state |A^), of the A-th excited state with A^ + 1 
particles |A^ -|- I, A) and of the I-th excited state with 
A^ — 1 particles jA^ — 1,/), respectively. The Lehmann 
representation of the one-particle Green's function is ob- 
tained by Fourier transform 



y^ /,(xi)/;(x2) 

^ oj-£i-iO+' 



rif \ Y^ /a(xi)/^(x2) 
G xi,X2;w = > ^ 

(3) 
where /a(x) = (7V|V'(x)|A^ 4- 1,A) and //(x) = (Af - 
1, /|V'(x)|A'') are the Dyson orbitals, and Ea = ^Af+i,yi — 
iSjv and El = E^ — E^^ij are minus the electron affini- 
ties and ionization energies, respectively. 



B. Two-particle Green's function 

The time-ordered two-particle Green's function is de- 
fined as 

^^02(1, 2; 1', 2') = (Ar|f[^(l)*(2)*^(2')^^(l')]|^^)- 

(4) 
Depending on the time ordering, it describes the propa- 
gation of a pair of holes, of electrons or of a hole and an 
electron. In the case of optical absorption, one is only 
interested in the propagation of a hole-electron pair. 

Let X be the 4-point polarizability, 

X(l, 2; 1', 2') = *G2(1, 2; l', 2') - *G(1, l')G(2, 2'). (5) 



where i^ = ii -I- 0+ . The Lehmann representation of the 
response function explicitly gives the excitation energies 
as poles in oj, 

x(xi,X2;xi,X2;w) = 

{N\¥{-x.{)^{-Ki)\N,K){N,K\¥{-}d^)^{y.2)\N) 



E 



^ — {En.k — E 



N 



iO^ 



E 



{N\'^^(li'2)i'{x2)\N,K){N,K\i'^l<L[)il{xl)\N) 



^ + {En.k - E_ 



N 



iO^ 



(7) 



where \N,K) is the K-th excited state with A^ par- 
ticles of energy En,k- The ground state \N,Q) = 
\N) is excluded from the sum. It is also use- 
ful to define the independent-particle (IP) polar- 



izability X7P(1,2;1',2') 



■iG(l,2')G(2,l'). Its 



Lehmann representation is easily obtained by calculat- 
ing x/p(xi,X2;x'i,x^;t) = -iG(xi,x^; t)G(x2,x'i; -r) 
with equation ([2]) and taking the Fourier transform 



X/p(xi,X2;xi,X2;a;) 



^ /;(x^i)/a(xi)/1(x'2)//(x2) 

lA 



E 

lA 



uj^{8a~ Si) + iO+ 
/;(x^)/^(x2)/l(x'i)/,(xi) 



UJ + {£a- £i) 



iO+ 



(8) 



In practice, the one-particle and two-particle Green's 
function can be calculated with equations of motion. 



C. Dyson equation 

To make easier the connection with expressions in a 
finite spin-orbital basis, we systematically use 4-point in- 
dexes for all the two-electron quantities. The starting 
point is therefore a fully non-local time-dependent Hamil- 
tonian, 

H{ti)^ f dxidl'i'''{l)h{l,l')i'{l') 

dxid2dl'd2'¥{l)¥{2)v{l, 2; l', 2')^(l')^(2'), 

(9) 

where t>(l, 2; 1', 2') = z;ee(|ri - r2|)5(ti, ^2)^(1, l')'5(2, 2') 
is the spin-independent instantaneous Coulomb electron- 
electron interaction and ft.(l, 1') is the one-electron 
Hamiltonian which contains the electron kinetic opera- 
tor and the nuclei-electron interaction Vne, 



It describes the coupled motion of two particles minus 
the motion of the independent ones. When the times are 
appropriately ordered, the 4-point polarizability reduces 
to the linear response function 

x(xi,X2;x'^,X2;t) = x(xi,ii,X2, ^2; x'^, t]^,X2, 4) (6) 



hil, 1') = -S{1, 1')^ + 5{1, l')Ke(ri). (10) 

Using the equations of motion for the Heisenberg cre- 
ation and annihilation operators in the expression of the 
derivative of G with respect to time 17[ , one can obtain 



the following equation, 
i /"d3(5(l,3)^G(3,2) 



d3/i(l,3)G'(3,2) 



+ i 



I dMl'di'v[l, 3; 1', 3')G'2(1', 3'+; 2, 3++) = (5(1, 2), 

(11) 



where ++ stands for i^ +0+. A whole series of equations 
can be derived for the Green's functions, relating the 
one-particle Green's function to the two-particle Green's 
function, the two-particle one to the three-particle one, 
etc. But solving this set of equations is not wanted. 

To avoid this, one can use the Schwinger derivative 
technique. Introducing an external time-dependent po- 
tential t/(l,l') ~ [/(xi,x']^, ti)(5(ti, t']^), one can express 
the two-particle Green's function in terms of the one- 
particle one and of its derivative with respect to C/, eval- 
uated at t/ = 0, 



JG(1,2) 



-G2(1,4;2,3) + G(1,2)G(4,3). 



(5t/(3,4) 

Using this relation in equation pip , one can get 
d 



(12) 



/ 



d3 



i(5(l, 3)^-/1(1, 3) 



G(3,2) 
- / d3SH,,(l,3)G(3,2)=<5(l,2), 



(13) 



where YiHxc{^i 2) is the Hartree-exchange-correlationself- 
energy which takes into account all the two-particle ef- 
fects. It can be decomposed into a Hartree contribution 

Sh(1, 2) = -i I d3d3'v{l, 3; 2, 3')G(3'+, 3++), (14) 

and an exchange-correlation one 

S.c(l,2) = 

i d3dl'd3'd4v{l,3;l',3' 



^5U{3++,3'+r ^ ' ^ 



(15) 

One can define a Green's function Gh which shows no 
two-particle effects and therefore follows the equation of 
motion 



d'S 



^S{l,3)^-h{l,3) 
oti 



Gh{3,2)^d{l,2). (16) 



Using this relation in equation (|13p . one finally gets the 
Dyson equation for the one-particle Green's function, 

J d3 [GZ\1,3) - ^H.c{l,3)] G(3, 2) - <5(1, 2). (17) 

This equation is also often used under the forms 

G(l, 2) = G,,(l, 2) + /" d3dAGhil, 3)2^,^(3, 4)G(4, 2), 

(18) 



or 



G-i(l,2) = G^i(l,2)-E^,,(l,2). 



(19) 



D. Bethe-Salpeter Equation 

Starting from the Dyson equation (J19p . and taking the 
derivative with respect to G, one can get the so-called 
Bethe-Salpeter equation (see, e.g., Ref. [35[) 

X-'ih 2; 1', 2') = x7p(1- 2; l', 2') - EhU^ 2; 1', 2'), 

(20) 
or 



X(l,2;l',2') = x/p(l,2;l',2') 



J d3456x/p(l, 4; 1', 3)2^^,(3, 6; 4, 5)x(5, 2; 6, 2'), 



(21) 



where ^h^hxc is the Hartree-exchange-correlation Bethe- 
Salpeter kernel, defined as 



•=ff2;c(3, 6;4, 5) = i 



<5G(5,6) ■ 



(22) 



E. Hedin's equations 

We now have equations of motion for the one- and two- 
particle Green's functions. They depend on the Hartree- 
exchange-correlation self-energy. Only its Hartree part 
is known exactly. A practical way of calculating its 
exchange-correlation part is needed. Hedin proposed a 
scheme which yields to a set of coupled equations and 
allows in principle for the calculation of the exact self- 
energy [3g. This scheme can be seen as a perturbation 
theory in terms of the screened interaction W instead of 
the bare Coulomb interaction v. We show a generaliza- 
tion of this derivation for the case of a non-local potential. 

Let F(5,6) = U{5,6) -iJd3d3'v{5,3;6,3')G{3',3+) 
be the non-local classical potential. Using the chain rule 
in the exchange-correlation self-energy, we get: 

S^c(l,2) = -i d3dl'd3'dAd5d6v{l,3;l' ,3')G{1' ,4) 

6G-^{A,2) 6V{5,6) 
^ 5V{5,6) 6U{3++,3'+) 

= i I d3dl'd3'd4d5d6v{l, 3; l', 3')G(l', 4) 

xf(4,6;2,5)e-^(5,3';6,3+). 

(23) 

where the inverse dielectric function e~^ which screens 
the bare Coulomb interaction v and the irreducible vertex 
function F are defined by 



£-1(1,2; 3,4) 



'^^(1,3) 

6U{A,2) 



and f(l,2;3,4) 



6G-\1,3) 
' 5V{A,2) ■ 



(24) 



We can therefore define a dynamically screened potential 
W{1, 2; 1', 2') = / d3d3'e~\l, 3; l', 3'+)w(2, 3'; 2', 3) 



= / d3d3'e~i(l, 3; 1', 3'+)w(3', 2; 3, 2') 



(25) 

where the symmetry of the Coulomb interaction v has 
been used, and we get the expression of the exchange- 
correlation self-energy, 

Sxc(l,2) = 

i I dl'd3d3'd4:G{l', 4)f(4, 3', 2, 3)M^(3, 1; 3', !')• 

(26) 

We still need to express the dielectric function and the 
irreducible vertex function without the use of V and U. 
To achieve this, we define the irreducible polarizability 
X(l, 2; 3, 4) = -iSG{l, 3)/SV{4, 2), which, with the prop- 
erties of the inverse and the definition of the vertex cor- 
rection, can be rewritten as 

X(l, 2; 3, 4) = -i / d5d5'G{l, 5)G(5', 3)f(5, 2; 5', 4). 

(27) 
Using this relation, one can rewrite the dielectric function 
as 

e(l,2;3,4) = 

5(1, 4)<5(2, 3) - / d5d5'v{l, 5; 3, 5')x(5', 2; 5+, 4), 

(28) 

and the irreducible vertex correction as 

f(l,2;3,4) = 

<5(1, 4)6(2, 3)-ij d5d6^g2ik^x(5, 2; 6, 4). ^^^^ 

We now have a set of five coupled equations ([25]) to ((29|) 
to calculate the self-energy. In practice, this set of equa- 
tions is never solved exactly, approximations are made. 

F. Static GW approximation 

We discuss now the static GW approximation which 
is the most often used approximation in practice in the 
BSE approach. 

In the GW approximation, one takes r(l,2;3,4) = 
(5(1,4)5(2,3). This simplifies greatly Hedin's equations. 
The irreducible polarizability becomes x(l,2;3,4) = 
-iG(l,4)G(2,3) = x/p(l,2;3,4) and the exchange- 
correlation self-energy becomes 

E,e(l, 2) = i I dl'd3G{l', 3)W{3, 1; 2, 1'). (30) 



If the derivative of W with respect to G is further 
neglected, as usually done, the corresponding Bethe- 
Salpeter kernel is then 

Sff,,(l, 2; 1', 2') = v{l, 2; l', 2') - W{2, 1; l', 2'), (31) 

where W is obtained from equation (|25[) and e~^ 
with equation p8|l in which x is replaced by xip. 
The Coulomb interaction is instantaneous and the one- 
particle Green's functions depends only of the time dif- 
ference, therefore the time dependence of the screened 
interaction is 

I^(l,2;l',2') = W(xi,X2;x;,x'2;r)5(ti,t;)5(t2,t2), 

(32) 
where t = ti — ^2- If one considers the time dependence in 
W , the Fourier transform of the Bethe-Salpeter equation 
is not straightforward [S^]- We will only consider the 
usual approximation where the screened interaction is 
static, i.e., 

iy(l,2;l',2')=W^(xi,X2;xl,x'2)5(ti,i'i)5(t2,4)'5(ti,i2). 

(33) 
To summarize, the Fourier-space Bethe-Salpeter equa- 
tion in the static GW approximation writes 

X~^(xi,X2;x3,X4;w) = 

X7p(xi, X2; Xg, X4; Uj) - S^3.c(Xi, X2; X3, X4), 

(34) 

where the kernel Si:/a;c(xi,X2; X3,X4) = 

w(xi,X2;x3,X4) — W^(x2,Xi;x3,X4) contains the static 
screened interaction W calculated from 

I^(xi,X2;Xi,X2) = 

/ (ix3dx3e~^(xi,X3;xi,X3)v(x3,X2;x3,X2) 



(35) 



and 



e(xi,X2;x3,X4) = (5(xi,X4)(5(x2,X3) 

- / dx5dx5v(xi,X5;x3,X5)x/p(x5,X2;x,5,X4;u; = 0). 

(36) 

We will refer to the approach of equations (|34p - ([5^ as the 
BSE-GW method. The one-particle Green's function G 
in Xip ~ —iGG is not yet specified. Different choices can 
be made. The simplest option is to use a non-interacting 
Green's function Go from a Hartree-Fock (HF) or Kohn- 
Sham (KS) calculation. In this case, xip ~ —iGoGo = 
Xo is just the non- interacting HF or KS response func- 
tion. In the condensed-matter physics literature, the 
usual recipe is to use xo in equation (|36p but an improved 
Xip in equation ((34|) from a GW calculation. In the 
case of H2 in a minimal basis, it is simple enough to use 
Xip constructed with the exact one-particle Green's func- 
tion G. Finally, we note that the dielectric function of 



equation ([36]) could be alternatively defined as including 
the HF exchange in addition to the Coulomb interaction, 
i.e. w(xi,X5;X3iX^) -^ ?;(xi,X5;x3,x[;)-i;(x5,Xi; X3,x^) 
(see, e.g., Ref. ISTl). which removes the "self-screening er- 
ror" for one-electron systems [3^ , but we do not explore 
this possibility here. 



III. EXPRESSIONS IN A FINITE ORBITAL 
BASIS 

A. Spin-orbital basis 

In order to solve the Bethe-Salpeter equation for fi- 
nite systems, all the equations are projected onto an or- 
thonormal spin-orbital basis {4>p}. As the equations are 
4-point equations relating two-particle quantities, they 
are in fact projected onto the basis of products of two 
spin orbitals. Each matrix element is thus indexed by 
two double indices. 

We consider the simplest case for which xip = Xo- 
The Lehmann representation of xo is 

Xo(xi,X2;xi,X2;w) = 

^0*(x'i)0,(xi)0:(x'2)0,(x2) 






{ea~e^) + iO+ 



(37) 



Ll!+ {Sa- £i) - ^0+ 



where c/ji is the i-th occupied spin-orbital of energy Si and 
(j)a is the a-th virtual spin-orbital of energy Sa- One can 
notice that xo is expanded only on occupied-virtual (ov) 
and virtual-occupied (vo) products of spin-orbitals. The 
matrix elements of xo are given by 

[Xo{^)]pq^rs = / rfxidx'irfx2dx^0p(x'i)0*(xi) 

X xo(xi,X2;x'i,X2;w)0*(x2)0s(x2). 

The matrix representation of its inverse, in the (ov,vo) 
subspace, is 



Xo"'H = - 



Ae 
Ae 



1 
-1 



(39) 



where Asiajb = A£ai,bj = {sa - £i)SijSab, where i,j re- 
fer to occupied spin-orbitals and a, b to virtual orbitals. 
The dimension of the matrix is thus 2MoMy x 2MoM„ 
where Mo and My are the numbers of occupied and vir- 
tual spin orbitals, respectively. To build the matrix x~^j 
one then needs to construct the matrix elements of the 
Bethe-Salpeter kernel '^hxc which are given by 



(5 



Hxc)pq,rs 



^pq,rs 



w, 



pr.qs 1 



(40) 



where u„ 



= (grips) are the usual two-electron inte- 
grals, and the matrix elements of W can be obtained 



from equation ((35 



pq,rs — I dxidx[dX2dx.2<l>p{x[)(j)*{xi) 

X Ty(xi,X2;x;,x'2)<?!);^(x2)<?!)s(x^) 

(ixi(ix']^(ix2(ix2(ix3(ix3(/)p (x']^ )(/)* (xi ) 

X e"^(xi,X3;xi,X3)u(x3,X2;x3,X2)^*(x2)(/)s(x2). 

(41) 



To decouple the common coordinates in e~^ and v, 
one can introduce two delta functions (5(x3,X4) and 
(5(x3,X4) and use the closure relations (5(x3,X4) = 
Etrt{^3)M^A) and <5(x^,x^) = EuM^3)K{0- By 
doing so, the matrix elements of v and e ^ appear explic- 
itly and we get 



yypq_rs — 7 ^ ^pq,tu'^iu,rs- 
tu 

Similarly, for the dielectric function, we have 



(42) 



'^pq,rs 



Op'pOqg 



Vpq,rs [Xo(w = 0)]^ 



(43) 



where the last equality comes from the fact that xo ha-s 
only diagonal elements. It can be seen that the static 
screened interaction consists of an infinite-order pertur- 
bation expansion in the Coulomb interaction, namely us- 
ing matrix notations, 

W = e-^ -v 

= v + v-xn{^ = ^)-v (44) 

+ V ■ xo{^ ^0)-v ■ xo(w ^0)-v + ..., 

the first term in this expansion corresponding to time- 
dependent Hartree-Fock (TDHF). The matrix represen- 
tation of the inverse of the interacting response function, 
in the (ov,vo) subspace, is then 



A B 

B* A* 



1 
-1 



with the matrices 



(45) 



(46a) 



(46b) 



The block structure of equation (j45]) is a consequence of 
the symmetry of the Coulomb interaction, Vqp,sr = v*g,^g, 
and of the static screened interaction, Wqs,pr = Wp^.^^. 
Moreover, the matrix A is Hermitian (because Viajb = 
^Aia ^'^i^ Wij^ab = 1^7*i6a) ^"^^ ^he matrix B is sym- 
metric (because Via^b] = Vjb.ai and Wib,aj = Wja.u)- 
The excitation energies a;„ are thus found by solving the 



usual linear-response pseudo-Hermitian eigenvalue equa- 
tion, just as in TDDFT, 



A B 

B* A* 






1 
-1 



Y.„ 



(47) 



whose solutions come in pairs: excitation energies aj„ 
with eigenvectors (Jf„,l^j), and de-excitation energies 
— aj„ with eigenvectors (1^*,X*). For real spin-orbitals 
and \i A^B and A — B are positive definite, the eigen- 
values are guaranteed to be real numbers and the pseudo- 
Hermitian eigenvalue equation (|47p can be transformed 
into a half-size symmetric eigenvalue equation [3|. 

If instead of starting from xoi one starts from x/p = 
—iGG with the exact one-particle Green's function G, the 
equations get more complicated since the matrix repre- 
sentation of xip is generally not diagonal and not only 
has contributions in the (ov,vo) subspace of spin-orbital 
products but also in the occupied-occupied (oo) and 
virtual-virtual (vv) subspace of spin-orbital products. 
The dimension of the matrices thus becomes M^ x Af ^ 
where M is the total number of (occupied and virtual) 
spin orbitals. In this case, the number of solutions of 
the response equations is generally higher than the num- 
ber of single excitations, and in particular double ex- 
citations might be obtained even without a frequency- 
dependent kernel. Spurious excitations can also be found. 
This is similar to what happens in linear-response TD- 
DMFT [IMl]. We will show this later in the case of H2 
in a minimal basis. 



B. Spin adaptation 



tained by 



The spin-adapted screened interaction is ob- 



yy-pq,rs — / ^ ^pq.tu'^tu.rs 



(50) 



where t and u refer to spatial orbitals, and the singlet 
dielectric function ^e^q^rs = tpU\,r\s\ + tpU\,risi is given 
by 

^Epg.rs = 5pr5qs - 2Vpq^rs [Xo(w = 0)]rs,rs ' (^1) 

The bottom line is that the linear-response eigenvalue 
equation (|T7l) fully decouples into a singlet eigenvalue 
equation 



lA iJB \ / iX„ \ 1 f 1 \ / iX„ 



IB* lA* ; I iY„ i - ^M -1 M 'Y„ 



^ia.jh ^'^ha.bj ^^ib.aji 

and a triplet eigenvalue equation 

1 



3A 3B 
3 TD* 3 4* 



Yn 



3v \ - ^^ \ Q _1 



with the matrices 

-^ia.jb — ^^ia.jb ^ij.ab-) 



(52) 



with the matrices 

^Aia^jb = Aeta,jb + '2Via,jb - Wij^ab, (53a) 



(53b) 



3y 



(54) 



(55a) 



We give now the expressions for spin-restricted closed- 
shell calculations. For four fixed spatial orbitals referred 
to asp, q, r, s, the Bethe-Salpeter kernel has the following 
spin structure 



^p-lq-l.r^s^ ^p^q]^^r],s]^ 



























"pigt- 




"pigt, 





(48) 



which can be brought to a diagonal form after rotation 
(see, e.g., Refs. Issl. l39l 1401) 



1^ 

' — 'pq,rs 
3^ 








-'pq.rs 






377 






■^pq,rs 





(49) 



^7^ 



•pq,rs 



2n 



with a spin-singlet term ^S 

three degenerate spin-triplet terms ''2. 



pq,rs 



Wpr,qs and 



'pq^rs 



-W, 



pr,qs • 



It has been used that the Coulomb interaction v and the 
screened interaction W are spin independent: Vpq^rs = 



^pT?Ti''tsT 



W, 



Vp-\q^,risi 



ptgt.'-tsT 



Vplql,r^s-t 
Wp^q^^risi 



Vpiqi,risi and 

Wplql,r^s^ = 



Bia.ib = -W, 



ib.aj • 



(55b) 



IV. EXAMPLE OF Ha IN A MINIMAL BASIS 

As a pedagogical example, we apply the BSE-GW 
method to the calculation of the excitation energies of 
H2 in a minimal basis consisting of two Slater basis func- 
tions, Lpa and ipb, centered on each hydrogen atom and 
with the same exponent C = 1. This is a closed-shell 
molecule, therefore all the calculations are done with spin 
adaptation in a spatial orbital basis. The molecular or- 
bitals are -01 = (Va + Vb)l\/'2.{1 + Sab) (symmetry (jg) 
and -02 = (Va - V3(,)/a/2(1 - Sab) (symmetry cr„) where 
Sab is the overlap between ipa and ipb- The matrix rep- 
resentations of all two-electron quantities in the space of 
spatial-orbital products are of the following form 



/ -Pii.ii -F'11,22 
-P22.11 ^22,22 



-pL2ai -Pl2,22 

\ -P2iai ^21,22 



^11,12 Al,21 
^22,12 ^22,21 



\ 



Pl2S2 -Pl2,21 
^21,12 -P21,21 / 



(56) 



and we refer to the upper left block as the (oo,vv) block, 
and to the bottom right block as the (ov,vo) block. All 
the values of the integrals as a function of the internu- 
clear distance R can be found in Ref. \4^ Note that, in 
the condensed-matter physics literature, a simplified ver- 
sion of H2 in a minimal basis with only on-site Coulomb 
interaction is often used under the name "half-filled two- 
site Hubbard model" (see, e.g., Refs.miH) ^. 



BSE-GW method using the non-interacting 
Green's function 



The simplest approximation in the BSE-GW method 
is to start from the non-interacting HF Green's function 
Go, leading to the non- interacting HF linear response 
function xip = —iGoGo = xo whose matrix representa- 
tion reads 



/O 




Xo(w) 








\ 



T" 



uj- Ae 





-1 



(57) 



Ae / 



where Ae = £2 — £1 is the difference between the energies 
of the molecular orbitals V'2 and i/'i ■ The non- interacting 
linear response function has non-vanishing matrix ele- 
ments only in the (ov,vo) block, but it will be necessary 
to consider the other blocks as well for the screened in- 
teraction W. The matrix of the Coulomb interaction is 



(58) 



where Jpq = {p(l\p<l) and Kpq = {pq\qp) are the usual 
Coulomb and exchange two-electron integrals over the 
molecular orbitals V'l and ip2 ■ The off-diagonal blocks of 
V are zero by symmetry for H2 in a minimal basis, but 
this is not the case in general. By matrix product and 
inversion, we get the static singlet dielectric matrix 



/Jn 


J12 








J12 


J22 














K12 


K12 


V 





K12 


K12, 



/I 
1 








\ 



Ae 






(59) 



which, in this case, is block diagonal with the (00, vv) 
block being the identity. By using its inverse, we finally 



get the static screened interaction matrix 



/Jll J12 
J12 J22 



W 




















"IU2" 



- 4ii:i2/Ae 
K12 



- 4ii:i2/Ae 
K12 



1 + 4Xi2/Ae 1 + 4Ki2/Ae ) 



(60) 

which is block diagonal and the (00, vv) block is just the 
bare Coulomb interaction in the case of H2 in a mini- 
mal basis, but this is not generally true. We have then 
everything to construct the ^A and ^B matrices of equa- 
tion (|52|) for singlet excitations, which in the present case 
arc just one-dimensional 



M = Ae + 2Xi2- J12, 



^B = 2A'i2 



K, 



l + 4Xi2/Ae' 



(61a) 



(61b) 



and the ^ A and ^B matrices of equation ([54]) for triplet 
excitations 



U = Ae- J12, 



'B 



K 



12 



1 + 4A'i2/A£ 



(62a) 



(62b) 



Solving then the response equations by the standard 
Casida approach |3| , we get the singlet excitation energy 



Ae + AK12 - Ji: 



K12 



X Ae - J 



K 



l + 4A'i2/Ae 
1/2 



12 



12 



1 + AKi2/Ae 
and the triplet excitation energy 



(63) 



Ae - J12 



Ae - J12 



A"l2 



1 + 4A:i2/Ae 



1 + 4/vi2/Ae 



1/2 



(64) 



Note that, for this simple system, the A terms have the 
usual TDHF or configuration interaction singles (CIS) 
forms, and the screening has an effect only on the B 
terms, decreasing the exchange integral ^12 by a fac- 
tor of 1 -I- 4:Ki2jAe. Therefore, in the Tamm-Dancoff 
approximation |44| . which consists in neglecting B, the 
effect of screening would be lost and the method would 
be equivalent to CIS. It is interesting to analyze the effect 
of the screening as a function of the internuclear distance 
R. For small R, the orbital energy difference Ae is much 
greater than the exchange integral i^i2, so the screening 
factor 1 -|- 4i^i2/Ae is close to 1 and TDHF excitation 
energies are recovered. For large R (dissociation limit). 
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FIG. 1: Excitation energies of the singlet ^I]+ (left) and the triplet '^E+ (right) states of H2 in a minimal basis as a 
function of the internuclear distance R calculated by FCI, TDHF, and BSE-GW with the non-interacting HF 

Green's function Gq. 



Ae goes to zero, so the screening factor diverges and the 
term A'i2/(1 + 4:Ki2/Ae) vanishes. 

The excitation energies from the ground state ^E+ to 
the first singlet ^E+ and triplet '^E+ excited states are 
plotted as a function of R in Figure [TJ The reference 
curves are from a full configuration-interaction (FCI) cal- 
culation giving the exact excitation energies in this ba- 
sis. In a minimal basis, the singlet ^S^ excited state 
is constrained to dissociate into the ionic configuration 
H~ ... H^, and so in the dissociation limit i? — > 00 
the exact singlet excitation energy goes to a constant, 
I — A Ki 0.625 hartrce where / and A are the ioniza- 
tion energy and electron affinity of the hydrogen atom. 
The triplet ^E+ dissociates into the neutral configura- 
tion H'... H', as does the ground state, and so the 
exact triplet excitation energy goes to zero in the dis- 
sociation limit. TDHF gives accurate excitation ener- 
gies for small R, but gives qualitatively wrong curves in 
the dissociation limit. For the singlet state, the TDHF 
excitation energy goes to zero, a wrong behavior in- 
herited from the vanishing Ae in this limit. For the 
triplet state, the TDHF response equation suffers from 
a triplet instability for R > 4 bohr and the excitation en- 
ergy becomes imaginary. It is known that TDDFT with 
standard density-functional approximations gives simi- 
larly incorrect energy curves [3, |42, I45l - l47| . The BSE- 
GW method using the non-interacting HF Green's func- 
tion Go gives accurate excitation energies at small _R, 
but fails in the dissociation limit. The singlet excita- 
tion energy becomes imaginary for R > 4.9 bohr. In- 
deed, in the dissociation limit, Ae goes to zero and equa- 
tion (|63p leads to a negative term under the square root: 
^oj -^ ^(4/ii2- Ji2)(-Ji2). Similarly, the BSE-GW 
triplet excitation energy is imaginary between R — 4.0 
and R = 4.9 bohr, and incorrectly tends to a non-zero 
value in the dissociation limit. 

The BSE-GW method using the non-interacting HF 
Green's function Go thus badly fails for H2 in the dis- 



sociation limit. As this method is based on a single- 
determinant reference, this should not come as a surprise. 
However, the BSE approach also allows one to start from 
an interacting Green's function G taking into account the 
multiconfigurational character of stretched H2 . We will 
now test this alternative approach. 



BSE-GW method using the exact Green's 
function 



1. Independent-particle response function 

We apply the BSE-GW equations (I34l)-(l36| with the 
independent-particle response function xip ^ —iGG 
constructed from the exact one-particle Green's func- 
tion G, and which can be calculated by the Lehmann 
formula JS]) using the iV-electron ground state and the 
{N ± l)-electron states. The states to consider for H2 in 
a minimal basis are given in Figure[5] The ground state is 
composed of two Slater determinants, its energy is -Bat = 
2ei — Jii + Ec where Ec = A- -yA^^fl?^ is the corre- 
lation energy with 2 A = 2Ae + Ju + J22 - 4 J12 + 2/-i'i2. 
The coefficients of the determinants are determined by 
C2 = ciKi2/{A + ^Kl^ + A2) and cl+cl^ 1. The en- 
ergies of the two [N + l)-clcctron states are: i?jv+i.i = 
2ei-l-e2- Jii and E^+ia = 2e2-|-ei - Ju -I- J22 -2Ji2 + 
A'i2. The energies of the two {N — l)-clectron states are: 
En-1,1 = ei - Jii and Em-i,2 = £2 - 2Ji2 + K12. We 
thus obtain four poles for the exact one-particle Green's 
function. Two of them correspond to minus the electron 
affinities, 
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FIG. 2: A^-electron ground state, and {N ± l)-clcctron states for H2 in a minimal basis. 



and the other two correspond to minus the ionization 
energies, 



£1 — En — En-1,1 = El + Ec 



(66a) 



£[=E, 



N 



En -1,2 = 2ei — £2 — 'hi + 2Ji2 — K12 + Ec- 

(66b) 
In condensed-matter physics, £1 and £2 are associated 
with "quasi-particle" peaks of photoelectron spectra, 
while £[ and £'2 are associated with "satellites". The 
Dyson orbitals are also easily calculated, and we finally 
arrive at the matrix representation of xip in the basis of 
the products of spatial orbitals 



X/p(w) 



( xiPAi(y>) 

XJP,22(^) 



V 



\ 



with the matrix elements 



XIPS2{^) 

X/P.2l(w) / 

(67) 



X/P,ii(^) 



X/P22(w) 



X/P,12('^) 



„2 2 
C]_C2 



^2^2 



uj-[£'2-£i) uj+{£'2-£iy 



„2 2 
C1C2 



„2 2 
C]_C2 



W- (52-^0 CJ + (52-f()' 



(68a) 



(68b) 



' f / \ ' 



oj-{£2-£i) uj + {£'2-£[ 



X/P,2l(w) = 



CJ-(5^-£:;) L0+{£2-£l)' 



(68d) 



Therefore, whereas Xo('^) ha-s only one positive pole, 
Xip{ijj) has four distinct positive poles (and four symmet- 
ric negative poles). These poles are plotted in Figure [3l 
The lowest one, £2 — £1, called fundamental gap in the 
condensed-matter physics literature, can be considered 
as an approximation to a neutral single excitation energy 
since in the limit of non-interacting particles it equals the 
difference of the orbital eigenvalues Ae = £2 — ei. The 
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FIG. 3: Positive poles of the independent-particle linear 

response function in function of the internuclear 

distance R 



two intermediate poles, £2 — £1 and £2 ^ ^1, can be inter- 
preted as approximations to a double excitation energy 
since they reduce to 2A£ in the limit of non-interacting 
particles. Surprisingly, the highest pole, £3 ~^ii reduces 
to 3A£ in this limit and it is thus tempting to associate it 
with a triple excitation even though the system contains 
only two electrons! In the dissociation limit R — >■ 00, the 
four poles tends to the same value, i.e. / — A « 0.625 
hartree which is also minus twice the correlation energy 
—2Ec, showing that the non- vanishing fundamental gap 
in this limit is a correlation effect. Note that it has been 
shown [38| that the non-self-consistent GW approxima- 
tion (Go Wo) to the one-particle Green's function gives a 
fundamental gap which is too small by a factor of 2 in 
the dissociation limit, so we do not consider this approx- 
imation here. 



2. Excitation energies 

Having calculated the independent-particle response 
function, the next steps of the BSE-GW calculation of the 
excitation energies proceed similarly as in Section llV Al 
even though the expressions get more complicated. From 
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FIG. 4: Excitation energy of the singlet ^I]+ state of H2 

in a minimal basis as a function of the internuclear 

distance R calculated by FCI and BSE-GW with the 

exact Green's function. The lowest pole of x/p(a;), the 

fundamental gap £2 ~ £1, is also plotted for comparison. 
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FIG. 5: Excitation energy of the second singlet ^E^ 

state of II2 in a minimal basis as a function of the 

internuclear distance R calculated by FCI and BSE-GW 

with the exact Green's function. The poles £2 ~ ^1 ^^^ 

£2 ~ £'i of xipi^) are also plotted for comparison. 



the matrix xip{^ ~ 0) and the Coulomb interaction 
matrix (|58p . we calculate the singlet dielectric matrix 
which is still block diagonal but the upper left block is 
no longer the identity matrix. We calculate then the 
static screened interaction matrix which is still block di- 
agonal but the elements of its upper left block are now 
also affected by screening. We can then construct the cor- 
responding singlet and triplet Bethe-Salpeter kernel ^H 
and '^H. The response eigenvalue equations (|52|) and ([54]) 
are no longer applicable, so the singlet excitation energies 
are found by searching the values of w giving vanishing 
eigenvalues of the inverse singlet linear response matrix 
^x{^)~^ — Xip{^)~^ ~ ^'^j and the triplet excitation en- 
ergies are found by searching the values of lo giving van- 
ishing eigenvalues of the inverse triplet linear response 
matrix '^x{<^)~^ = Xip{'^)~^ ~ '^'^- For II2 in a minimal 
basis, ^%(w)~^ and ^xi'^)^^ arc 4x4 matrices which are 
block diagonal, the (00, vv) block being uncoupled to the 
(ov,vo) block. For both the singlet and triplet cases, the 
four positive poles of xip{^) transform into four exci- 
tation energies (plus four symmetric de-excitation ener- 
gies). 

Among the two positive excitation energies coming 
from the (ov,vo) block of the matrix ^x(cj)~^, the low- 
est one is identified with the first singlet ^Sj excitation 
energy, which is called the optical gap. It is plotted in 
Figure |4] and compared with the reference FCI excitation 
energy and also with the fundamental gap £2 — ^1 to high- 
light the eflfect of the Bethe-Salpeter kernel. At small 
internuclear distance, R < 3 bohr, the Bethe-Salpeter 
kernel brings the BSE-GW curve is very close to the FCI 
curve. For large R, the BSE-GW excitation energy fol- 
lows the curve of the fundamental gap, which slightly 
overestimates the excitation energy at R ~ 10 bohr but 
eventually goes to the correct limit I — A when i? — > 00. 
Thus, contrary to the BSE-GW method using the non- 



interacting Green's function, the obtained excitation en- 
ergy curve has now a correct shape. This relies on the 
fundamental gap being a good starting approximation to 
the optical gap. As regards the second excitation energy 
coming from the (ov,vo) block of the matrix ^x(a;)~^ 
which is connected to highest pole £2 — £[ of x/p(^), it 
is a spurious excitation due to the approximate Bethe- 
Salpeter kernel used. 

The lowest positive excitation energy coming from the 
(00, vv) block of the matrix ^x{^)~^ is identified with the 
second singlet 



^E^ excited state which has a double exci- 
tation character. It is plotted in Figure [5] and compared 
with the FCI excitation energy for this state and with the 
poles £2 — £1 and £2 — £[ of xip {'^)- It is noteworthy that 
the BSE-GW method starting from xip{'^) instead of 
Xo('^) but using a frequency-independent kernel docs de- 
scribe this double-excitation state with an overall correct 
shape for the energy curve. However, the BSE-GW exci- 
tation energy is almost identical to the two poles £2 — £1 
and £2 — £[■ The Bethe-Salpeter kernel in the static GW 
approximation thus brings virtually no improvement for 
this state over the starting poles of x/p(w). The (00, vv) 
block of the matrix ^x(w)^^ also gives a second higher 
excitation energy that is spurious. 

We finally consider the triplet excited state '^E^ . The 
lowest positive excitation energy coming from the (ov,vo) 
block of the matrix ^x(aj)~^ should be identified with this 
state. It is plotted in Figure[n]and compared with the FCI 
excitation energy for this state and with the fundamental 
gap £2 — ^1- For small internuclear distances, i? < 3 bohr, 
the BSE-GW method gives an accurate excitation energy, 
but for larger R, instead of going to zero, the BSE-GW 
excitation energy follows the fundamental gap until the 
excitation energy becomes imaginary for R > 6.5 bohr. 
The problem is that the poles of xipi^) are the same for 
both the singlet and triplet cases, and the fundamental 
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FIG. 6: Excitation energy of the triplet ■^1]+ state of H2 

in a minimal basis as a function of the internuclear 

distance R calculated by FCI and BSE-GW with the 

exact Green's function. The lowest pole of x/p(a;), the 

fundamental gap £2 — £1, is also plotted for comparison. 



gap £2 — £1 is not a good starting approximation to the 
triplet excitation energy in the dissociation limit. The 
Bethe-Salpeter kernel in the static GW approximation is 
not able of compensating for this bad starting point. In 
addition to this excitation energy, the BSE-GW method 
gives three other spurious triplet excitation energies. 



CONCLUSION 
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ants give accurate excitation energies to the first singlet 
^E+ and triplet "^E^ excited states. In the dissocia- 
tion limit, however, the two variants differ. The first 
variant, starting from the non- interacting one-particle 
Green's function, badly fails in this limit for both the 
singlet and triplet states, giving imaginary excitation en- 
ergies. The second variant, starting from the exact one- 
particle Green's function, gives a qualitatively correct en- 
ergy curve for the singlet ^I]+ excited state up to the 
dissociation limit. This relies on the fact that the funda- 
mental gap (given by the one-particle Green's function) is 
a good starting approximation to the first singlet excita- 
tion energy. However, the same variant gives an incorrect 
energy curve for the triplet ^E+ excited state in the dis- 
sociation limit. In this case, the fundamental gap is a 
bad starting approximation to the first triplet excitation 
energy. 

The second BSE variant using the exact one-particle 
Green's function gives more excitation energies that the 
first BSE variant. Most of them arc spurious excita- 
tions due to the approximate Bethe-Salpeter kernel used. 
However, one of them can be identified with the excita- 
tion energy to the singlet ^E+ excited state which has a 
double excitation character. It is remarkable that such 
a double excitation can be described at all without us- 
ing a frequency-dependent kernel. However, the Bethe- 
Salpeter kernel in the static GW approximation is in- 
sufficient to describe accurately the energy curve of this 
state, even around the equilibrium distance. 



We have applied the BSE approach in the static GW 
approximation for the calculation of the excitation ener- 
gies on the toy model of H2 in a minimal basis. We have 
tested two variants for the starting one-particle Green's 
function: the non-interacting HE one and the exact one. 
Around the equilibrium internuclear distance, both vari- 
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